<<BACK

(Remi/David)

library(sf)
library(marmap)
library(raster)
library(tidyverse)
library(ggplot2)

IMAGINE!

Let’s imagine that you are interested in a species in a given area and wish to know as much as possible about it. But, you can’t go out in the field because funding is running short. What you do have, however, is a certain knowledge of the tools that are available to you and you are wondering how far they can take you. Let’s find out.

Defining species and area of interest

    # Let's select Atlantic cod as our species of interest
    sp <- 'Gadus morhua'

    # And the gulf of St. Lawrence as our study area
    #define the corners
    latmax <- 52.01312
    latmin <- 45.52399
    lonmax <- -55.73636
    lonmin <- -71.06333
    #
    # latmax <- 46.25
    # latmin <- 45.5
    # lonmax <- -61.25
    # lonmin <- -62.25

Describe your species

Fishbase

We’ll start with a description of the species. First, let’s see what fishbase has to offer. This online data repository, along with sealifebase, contains a lot of information on marine and aquatic species all over the world and is accessible through the package [rfishbase](https://github.com/ropensci/rfishbase)

    # Let's have a look at the information available our species ecology
    ecol <- rfishbase::ecology(sp)
    ecol <- cbind(colnames(ecol), t(ecol))
    rownames(ecol) <- NULL
    ecol <- ecol[ecol[,2] != 0, ] # remove 0
    ecol <- ecol[!is.na(ecol[,2]), ] # remove NAs

    knitr::kable(ecol, col.names = c('Descriptors', 'Attributes'))
Descriptors Attributes
autoctr 33
sciname Gadus morhua
StockCode 79
EcologyRefNo 1371
HabitatsRef 1371
Neritic -1
Intertidal -1
Oceanic -1
Estuaries -1
Herbivory2 mainly animals (troph. 2.8 and up)
HerbivoryRef 5743
FeedingType hunting macrofauna (predator)
FeedingTypeRef 5743
DietTroph 4.09
DietSeTroph 0.179
DietTLu 4.34
DietseTLu 0.72
DietRemark Troph of adults from 7 studies.
DietRef 26813
FoodTroph 4.29
FoodSeTroph 1
FoodRemark Trophic level estimated from a number of food items using a randomized resampling routine.
AddRems Opportunistic predator that forages mainly at dawn and dusk (Refs. 1371, 46189). Larvae feed mainly on zooplankton while juveniles prey predominantly on benthic crustaceans; adults feed mainly on zoobenthos and fish (Refs. 5743, 9604, 26813) including juvenile cod. Fish prey becomes more common in the diet with increasing body size (Refs. 1371, 89387). Adults may cover large distances during the feeding period (Ref. 89387).
Young cod are also pre yed upon by different fish species and octopus. Adult cod are prey items of top predators like sharks, rays, whales, dolphins, seals, and sea birds (Refs. 9023, 9581, 26954, 43651, 45735).
In the Baltic it grows up to 5 kg weight in 7-8 years; in the North Sea it reaches 8 kg in the same time span . Natural mortality for adults of both stocks is assumed to be around M=0.2, resulting in a mean adult life expectancy and mean duration of the reproductive phase of 5 years (Ref. 88171).
Parasites of the speci es include protozoans (trypanosome), myxosporidians, monogeneid, trematodes, cestodes, nematodes, acanthocephalan, hirudinid and copepods (Ref. 5951).
Schooling -1
SchoolingFrequency sometimes
SchoolingLifestage juveniles and adults
SchoolShoalRef 1371
AssociationsRemarks Generally considered a demersal fish although its habitat may become pelagic under certain hydrogrphic conditions, when feeding or spawning. There is some evidence that cod leave the bottom and school pelagically to spawn in preferred temperatures when bottom tempetatures are unsuitable. Gregarious during the day, forming compact schools that swim between 30-80 m above the bottom, and scatter at night (Ref. 1371). Schooling behavior may be adaptive for feeding. Reproductive behavior during spawning involves the circling of a female often by only one male per spawning bout (Ref. 86779).
SoftBottom -1
HardBottom -1
Rocky -1
SeaGrassBeds -1
Entered 2
Dateentered 1991-10-17T00:00:00.000Z
Modified 2374
Datemodified 2014-02-06T00:00:00.000Z
SpecCode 69

Taxize

We can also get it’s taxonomy rather easily. The package taxize allows you to extract and validate, among other things, the taxonomy of millions of species by accessing an important number of online databases accessible through their Application program interface (API).

  • Encyclopedia of Life (EOL) eol
  • axonomic Name Resolution Service tnrs
  • Integrated Taxonomic Information Service (ITIS) itis
  • Global Names Resolver (from EOL/GBIF) gnr
  • Global Names Index (from EOL/GBIF)
  • IUCN Red List
  • Tropicos (from Missouri Botanical Garden)
  • Theplantlist.org
  • Catalogue of Life
  • National Center for Biotechnology Information
  • CANADENSYS Vascan name search API
  • International Plant Names Index (IPNI)
  • World Register of Marine Species (WoRMS)
  • Barcode of Life Data Systems (BOLD)
  • Pan-European Species directories Infrastructure (PESI)
  • Mycobank
  • National Biodiversity Network (UK)
  • Index Fungorum
  • EU BON
  • Index of Names (ION)
  • Open Tree of Life (TOL)
  • World Register of Marine Species (WoRMS)
  • NatureServe
    # Start by exporting the taxonomy of our species
        taxo <- taxize::classification(sp, db = 'worms', verbose = FALSE)
        taxo[[1]]
## # A tibble: 10 × 3
##              name       rank     id
##             <chr>      <chr>  <int>
## 1        Animalia    Kingdom      2
## 2        Chordata     Phylum   1821
## 3      Vertebrata  Subphylum 146419
## 4   Gnathostomata Superclass   1828
## 5          Pisces Superclass  11676
## 6  Actinopterygii      Class  10194
## 7      Gadiformes      Order  10313
## 8         Gadidae     Family 125469
## 9           Gadus      Genus 125732
## 10   Gadus morhua    Species 126436
    # We can also extract the common or scientific names using sci2comm() & comm2sci(), respectively.
        taxize::sci2comm(sp, db = 'itis')
## 
## Retrieving data for taxon 'Gadus morhua'
## $`Gadus morhua`
## [1] "morue de l'Atlantique"    "bacalao del Atl<e1>ntico"
## [3] "cod"                      "rock cod"                
## [5] "morue franche"            "Atlantic cod"
    # Or find out whether there are other names under which the species is known
        taxize::synonyms(sp, db = 'itis')
## 
## Retrieving data for taxon 'Gadus morhua'
## $`Gadus morhua`
##   sub_tsn acc_tsn       message
## 1  164712  164712 no syns found
## 
## attr(,"class")
## [1] "synonyms"
## attr(,"db")
## [1] "itis"
    # Another really interesting feature is to extract all known species at a given taxonomic scale.
    # Let's take the genus level and see the first 20 species that are part of that genus on the itis database
        taxize::children(strsplit(sp, ' ')[[1]][1], db = 'itis')[[1]]
## 
## Retrieving data for taxon 'Gadus'
## # A tibble: 4 × 5
##   parentname parenttsn rankname           taxonname    tsn
## *      <chr>     <chr>    <chr>               <chr>  <chr>
## 1      Gadus    164710  Species Gadus macrocephalus 164711
## 2      Gadus    164710  Species        Gadus morhua 164712
## 3      Gadus    164710  Species          Gadus ogac 164717
## 4      Gadus    164710  Species Gadus chalcogrammus 934083

GloBI

How about extracting all known preys and predators of the species of interest? The Global Biotic Interactions web platform contains thousands of empirical binary interactions for multiple types of interactions, all over the world, and is easily accessible using the package rglobi.

    # There are multiple types of interactions available on GloBI
        rglobi::get_interaction_types()[,1:3]
##           interaction     source     target
## 1                eats   consumer       food
## 2             eatenBy       food   consumer
## 3             preysOn   predator       prey
## 4        preyedUponBy       prey   predator
## 5               kills     killer     victim
## 6            killedBy     victim     killer
## 7          parasiteOf   parasite       host
## 8         hasParasite       host   parasite
## 9              hostOf       host   symbiont
## 10            hasHost   symbiont       host
## 11         pollinates pollinator      plant
## 12       pollinatedBy      plant pollinator
## 13         pathogenOf   pathogen       host
## 14        hasPathogen       host   pathogen
## 15           vectorOf     vector   pathogen
## 16          hasVector   pathogen     vector
## 17  dispersalVectorOf     vector       seed
## 18 hasDispersalVector       seed     vector
## 19         symbiontOf     source     target
## 20   flowersVisitedBy      plant    visitor
## 21    visitsFlowersOf    visitor      plant
## 22      interactsWith     source     target
    # For now let's focus on predator-prey interactions
        prey <- rglobi::get_prey_of(sp)$target_taxon_name
        pred <- rglobi::get_predators_of(sp)$target_taxon_name
        prey
##  [1] "Neocalanus tonsus"             "Pseudocalanus elongatus"      
##  [3] "Arctica islandica"             "Gastrosaccus spinifer"        
##  [5] "Diastylis rathkei"             "Buccinum undatum"             
##  [7] "Corystes cassivelanus"         "Eledone cirrhosa"             
##  [9] "Gonatus fabricii"              "Bathypolypus arcticus"        
## [11] "Rossia moelleri"               "Bathypolypus bairdii"         
## [13] "Todarodes sagittatus"          "Cancer pagurus"               
## [15] "Rossia macrosoma"              "Pandalus borealis"            
## [17] "Pandalus montagui"             "Lithodes maja"                
## [19] "Hyas coarctatus"               "Crangon allmanni"             
## [21] "Galathea strigosa"             "Ophiopholis aculeata"         
## [23] "Natatolana borealis"           "Atelecyclus rotundatus"       
## [25] "Crangon crangon"               "Carcinus maenas"              
## [27] "Mya arenaria"                  "Pagurus bernhardus"           
## [29] "Ophiothrix fragilis"           "Psammechinus miliaris"        
## [31] "Trisopterus luscus"            "Callionymus lyra"             
## [33] "Actinopterygii"                "Gadiformes"                   
## [35] "Gobiidae"                      "Eutrigla gurnardus"           
## [37] "Trachurus trachurus"           "Hippoglossoides platessoides" 
## [39] "Scomber scombrus"              "Lepidorhombus whiffiagonis"   
## [41] "Trisopterus esmarkii"          "Sardina pilchardus"           
## [43] "Trisopterus minutus"           "Enchelyopus cimbrius"         
## [45] "Buglossidium luteum"           "Microchirus variegatus"       
## [47] "Merlangius merlangus"          "Actinopterygii"               
## [49] "Ammodytes tobianus"            "Clupea harengus"              
## [51] "Mallotus villosus"             "Lethenteron camtschaticum"    
## [53] "Lumpenus lampretaeformis"      "Sprattus sprattus"            
## [55] "Pomatoschistus minutus"        "Zoarces viviparus"            
## [57] "Gadus morhua"                  "Ammodytes marinus"            
## [59] "Buenia jeffreysii"             "Lophius piscatorius"          
## [61] "Sebastes viviparus"            "Microstomus kitt"             
## [63] "Agonus cataphractus"           "Solea solea"                  
## [65] "Merluccius merluccius"         "Micromesistius poutassou"     
## [67] "Salmo salar"                   "Myxine glutinosa"             
## [69] "Glyptocephalus cynoglossus"    "Scophthalmus rhombus"         
## [71] "Pleuronectes platessa"         "Pollachius virens"            
## [73] "Scophthalmus maximus"          "Zeugopterus punctatus"        
## [75] "Melanogrammus aeglefinus"      "Phrynorhombus norvegicus"     
## [77] "Squalus acanthias"             "Merluccius bilinearis"        
## [79] "Zoarces americanus"            "Pseudopleuronectes americanus"
## [81] "Scophthalmus aquosus"          "Ammodytes dubius"             
## [83] "Limanda limanda"               "Myoxocephalus scorpius"       
## [85] "Scyliorhinus canicula"
        pred
##  [1] "Larus"                        "Thalasseus sandvicensis"     
##  [3] "Phalacrocorax carbo"          "Phoca vitulina"              
##  [5] "Sprattus sprattus"            "Clupea harengus"             
##  [7] "Scomber scombrus"             "no name"                     
##  [9] "Gadus morhua"                 "Thunnus thynnus"             
## [11] "Petromyzon marinus"           "Pollachius virens"           
## [13] "Molva molva"                  "Squalus acanthias"           
## [15] "Hippoglossoides platessoides" "Pomatomus saltatrix"         
## [17] "Myxine glutinosa"             "Lophius piscatorius"         
## [19] "Eutrigla gurnardus"           "Hippoglossus hippoglossus"   
## [21] "Amblyraja radiata"            "Anarhichas lupus"            
## [23] "Reinhardtius hippoglossoides" "Sebastes"                    
## [25] "Sebastes mentella"            "Xiphias gladius"             
## [27] "Merlangius merlangus"

Make your search spatially explicit

Get data from OBIS

OBIS is the Ocean Biogeographic Information System and their vision is: “To be the most comprehensive gateway to the world’s ocean biodiversity and biogeographic data and information required to address pressing coastal and world ocean concerns.”

We can get access to their HUGE database through the excellent robis package by ropensci (unsurprisingly, they also have other great open science R packages, from which this whole section is highly inspired)

    sf::st_as_text(bb$geometry)
    OBIS <- robis::occurrence(scientificname = sp, geometry=st_as_text(bb$geometry), startdate = as.Date("2010-01-01"), enddate = as.Date("2017-01-01"), fields = c("species", "yearcollected","decimalLongitude", "decimalLatitude"))

    # remove duplicates
        OBIS <- unique(OBIS)

    write.csv(OBIS,'obis.csv',row.names = FALSE)

Visualize data

We can easily visualize the data…

OBIS <- read.csv("obis.csv")
plot(OBIS[,c("decimalLongitude", "decimalLatitude")])

… but the aesthetics could be much better!

Plot your study site

Define site map limits

Let’s start by defining a bounding box for your study site:

#create a matrix
bbmat <- cbind(
            c(lonmin,lonmax,lonmax,lonmin,lonmin),
            c(latmin,latmin,latmax,latmax,latmin)
        )
# make the matrix a 'simple features' polygon
bbsf <- st_polygon(list(bbmat))
# and let's give it information about the projection:
bbsfc <- st_sfc(bbsf,crs="+proj=longlat +datum=WGS84")
# finally, let's make it a simple features data.frame
bb <- st_sf(name="Study Site",geometry=bbsfc)
str(bb)
## Classes 'sf' and 'data.frame':   1 obs. of  2 variables:
##  $ name    : chr "Study Site"
##  $ geometry:sfc_POLYGON of length 1; first list element: List of 1
##   ..$ : num [1:5, 1:2] -71.1 -55.7 -55.7 -71.1 -71.1 ...
##   ..- attr(*, "class")= chr  "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA
##   ..- attr(*, "names")= chr "name"
# plot it!
ggplot(bb) + geom_sf()

Add a basemap!

There are many ways to add a basemap in R

marmap package

bathymetry <- getNOAA.bathy(lonmin,lonmax,latmin,latmax,resolution=1,keep=TRUE)
## File already exists ; loading 'marmap_coord_-71.06333;45.52399;-55.73636;52.01312_res_1.csv'
plot.bathy(bathymetry,image=T,drawlabels=TRUE)

blues <- colorRampPalette(c("darkblue", "lightblue1"))
greys <- colorRampPalette(c(grey(0.4),grey(0.99)))

plot.bathy(bathymetry,
           image = TRUE,
           land = TRUE,
           n=0,
           bpal = list(c(0, max(bathymetry), greys(100)),
                       c(min(bathymetry), 0, blues(100))))

# this is a 'bathy' object and you could plot it as is with the tools provided in marmap, but to keep things 'simple' I'm going to convert to sf. There is no direct method yet, so... don't judge me!
# convert from bathy to raster to 'sp' polygon to simple features
# bathypoly <- marmap::as.raster(bathymetry) %>%
    # rasterToPolygons() %>% st_as_sf
bathyras <- fortify.bathy(bathymetry)

bathyras$z[bathyras$z>0] <- NA

x=0.2
ggplot() +
    geom_raster(data=bathyras,aes(x=x,y=y,fill=z))+
    geom_sf(data=bb,fill="transparent",colour='yellow',size=2)+
    # geom_sf(data=Canada,fill='grey40')+
    coord_sf(xlim = c(lonmin-x,lonmax+x),
             ylim = c(latmin-x,latmax+x),
             expand = F)

Overlay your species occurrences over the base map

You can now overlay your species occurrences on your base map!

OBIS <- st_as_sf(OBIS,
                 coords = c("decimalLongitude", "decimalLatitude"),
                 crs="+proj=longlat +datum=WGS84",
                 remove=FALSE) %>%
    filter(!is.na(species))

    ggplot(OBIS) +
        geom_raster(data=bathyras,aes(x=x,y=y,fill=z))+
        geom_sf(data=bb,fill="transparent",colour='yellow',size=2)+
        # geom_sf(data=Canada,fill='grey40')+
        geom_point(data=OBIS,aes(x=decimalLongitude,y=decimalLatitude))+
        facet_wrap(~species)+
        coord_sf(xlim = c(lonmin-x,lonmax+x),
                 ylim = c(latmin-x,latmax+x),
                 expand = F)

Environmental variables

There are also multiple resources to access environmental data. As seen before, we can easily access bathymetry data using the marmap package. Other environmental data can also be accessed using other packages like raster for terrestrial climatic data and rnoaa

Extract sea surface temperature (SST)

    # Bear with us, for simplicity's sake we will only take one year of data for one month, knowing very well that this does not make sense ecologically!
    # Eventually use apply() to load multiple datasets and create rasterStack
        # An exemple with sea ice
            # urls <- seaiceeurls(mo='Apr', pole='N')[1:10]
            # out <- lapply(urls, seaice)
            # names(out) <- seq(1979,1988,1)
            # df <- ldply(out)
            # library('ggplot2')
            # ggplot(df, aes(long, lat, group=group)) +
            # geom_polygon(fill="steelblue") +
            # theme_ice() +
            # facet_wrap(~ .id)

        sst.nc <- rnoaa::ersst(2010, 8, overwrite = TRUE)
        sst <- raster(sst.nc$filename, varname = 'sst')
        sst@data@values <- values(sst) # make sure values are in the right data slot

        # never do that...
        sst@extent@xmin <- sst@extent@xmin - 361
        sst@extent@xmax <- sst@extent@xmax - 361

    # Transform rasterLayer as a data.frame
        sst.spdf <- as(sst, "SpatialPixelsDataFrame")
        sst.df <- as.data.frame(sst.spdf)
        colnames(sst.df) <- c('sst','x','y')

    # Plot newly acquired data!
        ggplot(sst.df, aes(x=x, y=y)) +
            geom_tile(aes(x=x,y=y,fill=sst)) +
            coord_equal()

Species distribution models

Now what we can do once we have species occurrence data and environmental variables are species distribution models (SDMs). There are multiple ways and packages to perform SDMs and this vignette provides an extensive overview of the different methods available. For this workshop, we will be using the package biomod2, which implements most of the methods presented in the vignette. There is also a thorough explanation of our to use biomod2 in another very interesting vignette

Formating data

The first thing to do is to format our data to be able to work with biomod2 functions, which work with BIOMOD.formated.data.PA class objects. This entails stacking all of our environmental data together to perform analyses.

    # Create raster stack with environmental variables
    # raster have to be of similar extent and resolution
        # Load bathy data from marmap as raster
            bathypoly <- marmap::as.raster(bathymetry)
            bathypoly@data@values[bathypoly@data@values > 0] <- NA

        # Resample sst to make it same extent as bathypoly
            # create empty raster with appropriate resolution and extent
                sst2 <- raster(resolution = res(bathypoly), ext = extent(bathypoly), crs=bathypoly@crs)
            # Beware, this may not be good practice because we are changing the resolution of the data to a finer resolution!
                sst2 <- resample(sst, sst2, method='bilinear') # extract data from sst

        # Create raster stack
            envCov <- raster::stack(bathypoly, sst2)

        plot(envCov)

You can then create your biomod2 data class.

    OBIS <- read.csv("obis.csv")
    SDMdata <- biomod2::BIOMOD_FormatingData(resp.var = rep(1,nrow(OBIS)),
                                            expl.var = envCov,
                                            resp.xy = OBIS[, c("decimalLongitude","decimalLatitude")],
                                            resp.name = sp,
                                            PA.nb.rep = 2)
## 
## -=-=-=-=-=-=-=-=-=-=-=-= Gadus morhua Data Formating -=-=-=-=-=-=-=-=-=-=-=-=
## 
##  Response variable name was converted into Gadus.morhua
##       ! No data has been set aside for modeling evaluation
##    > Pseudo Absences Selection checkings...
##    > random pseudo absences selection
##    > Pseudo absences are selected in explanatory variables
##          ! Some NAs have been automaticly removed from your data
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= Done -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    SDMdata
## 
## -=-=-=-=-=-=-=-=-=-=-=-= 'BIOMOD.formated.data.PA' -=-=-=-=-=-=-=-=-=-=-=-=
## 
## sp.name =  Gadus.morhua
## 
##   1187 presences,  0 true absences and  1929 undifined points in dataset
## 
## 
##   2 explanatory variables
## 
##      layer         Extended.reconstructed.sea.surface.temperature
##  Min.   :-527.00   Min.   :11.46                                 
##  1st Qu.:-215.25   1st Qu.:15.63                                 
##  Median : -95.50   Median :16.67                                 
##  Mean   :-144.91   Mean   :16.33                                 
##  3rd Qu.: -54.75   3rd Qu.:17.15                                 
##  Max.   :   0.00   Max.   :18.12                                 
## 
## 
##  2 Pseudo Absences dataset available ( PA1 PA2 ) with  974 
## absences in each (true abs + pseudo abs)
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    plot(SDMdata)

You can then use this object with the biomod2 functions are generate a distribution model for your species of interest.

    # Basic options for modelling
        SDMoption <- biomod2::BIOMOD_ModelingOptions()

    # SDM model
        SDM <- biomod2::BIOMOD_Modeling(data = SDMdata,
                                        models = "MAXENT.Tsuruoka",
                                        model.options = SDMoption)
## 
## 
## Loading required library...
## 
## Checking Models arguments...
## Warning in .Models.check.args(data, models, models.options, NbRunEval,
## DataSplit, : Models will run with 'defaults' parameters
## 
## Creating suitable Workdir...
## 
##          ! Weights where automaticly defined for Gadus.morhua_PA1 to rise a 0.5 prevalence !
##          ! Weights where automaticly defined for Gadus.morhua_PA2 to rise a 0.5 prevalence !
## 
## 
## -=-=-=-=-=-=-=-=-=-=-= Gadus.morhua Modeling Summary -=-=-=-=-=-=-=-=-=-=-=
## 
##  2  environmental variables ( layer Extended.reconstructed.sea.surface.temperature )
## Number of evaluation repetitions : 1
## Models selected : MAXENT.Tsuruoka 
## 
## Total number of model runs : 2 
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
## 
## 
## -=-=-=- Run :  Gadus.morhua_PA1 
## 
## 
## -=-=-=--=-=-=- Gadus.morhua_PA1_Full
    # Print model summary
        SDM
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-= BIOMOD.models.out -=-=-=-=-=-=-=-=-=-=-=-=-=-=
## 
## Modeling id : 1493514294
## 
## Species modeled : Gadus.morhua
## 
## Considered variables : layer 
## Extended.reconstructed.sea.surface.temperature
## 
## 
## Computed Models :  Gadus.morhua_PA1_Full_MAXENT.Tsuruoka 
## Gadus.morhua_PA2_Full_MAXENT.Tsuruoka
## 
## 
## Failed Models :  none
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=
    # print the dimnames of this object to understand the structure of the returned objects
        dimnames(SDM)
## NULL
    # get all models evaluation.
    # These give you an evaluation of the performance of your model.
        biomod2::get_evaluations(SDM)
## , , MAXENT.Tsuruoka, Full, PA1
## 
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.192  439.0      91.997      26.078
## TSS          0.183  490.0      84.162      34.086
## ROC          0.581  490.5      84.078      34.189
## 
## , , MAXENT.Tsuruoka, Full, PA2
## 
##       Testing.data Cutoff Sensitivity Specificity
## KAPPA        0.190  464.0      87.616      30.021
## TSS          0.181  464.0      87.616      30.021
## ROC          0.583  465.5      87.616      30.538
    # Finally, once the models are calibrated and evaluated, you can project your species spatial distribution...
        SDMproj <- biomod2::BIOMOD_Projection(modeling.output = SDM,
                                              new.env = envCov,
                                              proj.name = 'Awesome!',
                                              selected.models = 'all',
                                              binary.meth = 'TSS',
                                              compress = 'xz',
                                              clamping.mask = F,
                                              output.format = '.grd')
## 
## -=-=-=-=-=-=-=-=-=-=-=-=-= Do Models Projections -=-=-=-=-=-=-=-=-=-=-=-=-=
## 
##  > Building clamping mask
## 
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer')
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer')
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer')
## *** in setMethod('BinaryTransformation', signature(data='RasterLayer')
    # ... and visualize them!                
        plot(SDMproj)

Riveting! We starting with a basic idea and ended up with a species description, an idea of the distribution of certain environmental parameters in our area of interest (could be better, we know!) and the potential spatial distribution of our species!

…but should we stop there?

The reality is, there is so much more going on in our area of interest than the spatial distribution of a single species. Let’s see what other insights we an get from available open source resources.

Community level analyses

One of the usual challenges in ecology is to perform analyses at the community scale often due to the difficulty of sampling multiple species in a single system. The reality is that our species of interest is very likely to interact and be impacted by other species.

Extract all occurrence known occurrence data in area of interest

OBIS once again gives us the opportunity to access multiple species occurrences. Let’s go back to OBIS and extract every available data in our system.

    # Load all known species occurrences from our area of interest
        multiOBIS <- robis::occurrence(geometry=st_as_text(bb$geometry), startdate = as.Date("2010-01-01"), enddate = as.Date("2017-01-01"), fields = c("species", "yearcollected","decimalLongitude", "decimalLatitude"))
        multiOBIS <- multiOBIS[!is.na(multiOBIS[,'species']), ] # remove species == NA

    # Remove duplicates
        multiOBIS <- unique(multiOBIS)

    # Select species for which there is at least n observations
        n = 50
        spOK <- names(which(table(multiOBIS[,'species']) >= n)) # id of species with at least n observations
        multiOBIS <- multiOBIS[multiOBIS[, 'species'] %in% spOK, ] # Subset of multiOBIS data with species with at least n observations
        multiOBIS[,'species'] <- as.factor(multiOBIS[,'species'])
        multiOBIS[,'yearcollected'] <- as.factor(multiOBIS[,'yearcollected'])
        write.csv(multiOBIS,'multiOBIS.csv',row.names = FALSE)
    # Read occurrence data file
        multiOBIS <- read.csv("multiOBIS.csv")

    # Structure of object
        str(multiOBIS)
## 'data.frame':    37291 obs. of  4 variables:
##  $ species         : chr  "Themisto compressa" "Themisto compressa" "Themisto compressa" "Themisto compressa" ...
##  $ yearcollected   : int  2010 2010 2010 2010 2010 2010 2010 2010 2010 2010 ...
##  $ decimalLongitude: num  -68.6 -68.2 -64 -63.7 -60.2 ...
##  $ decimalLatitude : num  48.7 48.8 47.7 47.7 47 ...
    # Visualize species list
        spList <- unique(multiOBIS[,'species'])
        length(spList)
## [1] 167
        spList
##   [1] "Themisto compressa"               
##   [2] "Themisto abyssorum"               
##   [3] "Temora longicornis"               
##   [4] "Scina borealis"                   
##   [5] "Scolecithricella minor"           
##   [6] "Parasagitta elegans"              
##   [7] "Paraeuchaeta norvegica"           
##   [8] "Triconia borealis"                
##   [9] "Oithona similis"                  
##  [10] "Oithona atlantica"                
##  [11] "Metridia lucens"                  
##  [12] "Metridia longa"                   
##  [13] "Eukrohnia hamata"                 
##  [14] "Dimophyes arctica"                
##  [15] "Calanus hyperboreus"              
##  [16] "Calanus glacialis"                
##  [17] "Calanus finmarchicus"             
##  [18] "Boreomysis arctica"               
##  [19] "Aglantha digitale"                
##  [20] "Acartia (Acartiura) longiremis"   
##  [21] "Reinhardtius hippoglossoides"     
##  [22] "Pandalus borealis"                
##  [23] "Myoxocephalus scorpius"           
##  [24] "Merluccius bilinearis"            
##  [25] "Hippoglossus hippoglossus"        
##  [26] "Glyptocephalus cynoglossus"       
##  [27] "Enchelyopus cimbrius"             
##  [28] "Chionoecetes opilio"              
##  [29] "Menidia menidia"                  
##  [30] "Gadus morhua"                     
##  [31] "Fundulus heteroclitus"            
##  [32] "Hyas araneus"                     
##  [33] "Skeletonema costatum"             
##  [34] "Pseudo-nitzschia delicatissima"   
##  [35] "Homarus americanus"               
##  [36] "Carcinus maenas"                  
##  [37] "Gasterosteus aculeatus"           
##  [38] "Apeltes quadracus"                
##  [39] "Gasterosteus wheatlandi"          
##  [40] "Palaemonetes vulgaris"            
##  [41] "Pseudopleuronectes americanus"    
##  [42] "Osmerus mordax"                   
##  [43] "Fundulus diaphanus"               
##  [44] "Pungitius pungitius"              
##  [45] "Pleuronectes putnami"             
##  [46] "Cancer plebejus"                  
##  [47] "Microgadus tomcod"                
##  [48] "Syngnathus fuscus"                
##  [49] "Tautogolabrus adspersus"          
##  [50] "Myoxocephalus aenaeus"            
##  [51] "Scophthalmus aquosus"             
##  [52] "Morone saxatilis"                 
##  [53] "Urophycis tenuis"                 
##  [54] "Alosa pseudoharengus"             
##  [55] "Brisaster fragilis"               
##  [56] "Ctenodiscus crispatus"            
##  [57] "Mallotus villosus"                
##  [58] "Triglops murrayi"                 
##  [59] "Hyas coarctatus"                  
##  [60] "Myoxocephalus octodecemspinosus"  
##  [61] "Leptoclinus maculatus"            
##  [62] "Illex illecebrosus"               
##  [63] "Aspidophoroides monopterygius"    
##  [64] "Liparis gibbus"                   
##  [65] "Icelus spatula"                   
##  [66] "Artediellus atlanticus"           
##  [67] "Hippoglossoides platessoides"     
##  [68] "Lycodes lavalaei"                 
##  [69] "Limanda ferruginea"               
##  [70] "Gymnocanthus tricuspis"           
##  [71] "Artediellus uncinatus"            
##  [72] "Clupea harengus"                  
##  [73] "Lithodes maja"                    
##  [74] "Aspidophoroides olrikii"          
##  [75] "Amblyraja radiata"                
##  [76] "Ammodytes dubius"                 
##  [77] "Eumesogrammus praecisus"          
##  [78] "Chlamys islandica"                
##  [79] "Lycodes vahlii"                   
##  [80] "Hemitripterus americanus"         
##  [81] "Malacoraja senta"                 
##  [82] "Arctozenus risso"                 
##  [83] "Lumpenus lampretaeformis"         
##  [84] "Leptagonus decagonus"             
##  [85] "Eumicrotremus spinosus"           
##  [86] "Myxine glutinosa"                 
##  [87] "Gadus ogac"                       
##  [88] "Phycis chesteri"                  
##  [89] "Gymnelus viridis"                 
##  [90] "Scomber scombrus"                 
##  [91] "Cryptacanthodes maculatus"        
##  [92] "Nezumia bairdii"                  
##  [93] "Melanostigma atlanticum"          
##  [94] "Cyclopterus lumpus"               
##  [95] "Lophius americanus"               
##  [96] "Anarhichas lupus"                 
##  [97] "Centroscyllium fabricii"          
##  [98] "Careproctus reinhardti"           
##  [99] "Lycenchelys verrillii"            
## [100] "Phocoena phocoena"                
## [101] "Boreogadus saida"                 
## [102] "Argentina silus"                  
## [103] "Pandalus montagui"                
## [104] "Hippasteria phrygiana"            
## [105] "Strongylocentrotus droebachiensis"
## [106] "Psilaster andromeda"              
## [107] "Solaster endeca"                  
## [108] "Cucumaria frondosa"               
## [109] "Crossaster papposus"              
## [110] "Lebbeus polaris"                  
## [111] "Sebastes mentella"                
## [112] "Pontophilus norvegicus"           
## [113] "Oceanites oceanicus"              
## [114] "Morus bassanus"                   
## [115] "Uria aalge"                       
## [116] "Uria lomvia"                      
## [117] "Fulmarus glacialis"               
## [118] "Rissa tridactyla"                 
## [119] "Larus marinus"                    
## [120] "Alle alle"                        
## [121] "Larus argentatus"                 
## [122] "Alca torda"                       
## [123] "Fratercula arctica"               
## [124] "Oceanodroma leucorhoa"            
## [125] "Puffinus griseus"                 
## [126] "Puffinus gravis"                  
## [127] "Hemiselmis virescens"             
## [128] "Quadricilia rotundata"            
## [129] "Amphidinium sphenoides"           
## [130] "Ceratoneis closterium"            
## [131] "Lennoxia faveolata"               
## [132] "Chaetoceros convolutus"           
## [133] "Pyramimonas plurioculata"         
## [134] "Telonema subtile"                 
## [135] "Plagioselmis prolonga"            
## [136] "Leucocryptos marina"              
## [137] "Pseudopedinella pyriformis"       
## [138] "Heterocapsa rotundata"            
## [139] "Gyrodinium fusus"                 
## [140] "Leptocylindrus minimus"           
## [141] "Gyrodinium flagellare"            
## [142] "Protoperidinium bipes"            
## [143] "Lohmaniella oviformis"            
## [144] "Gymnodinium verruculosum"         
## [145] "Gyrodinium formosum"              
## [146] "Gymnodinium galeatum"             
## [147] "Torodinium robustum"              
## [148] "Teleaulax amphioxeia"             
## [149] "Chaetoceros similis"              
## [150] "Monosiga marina"                  
## [151] "Rhodomonas marina"                
## [152] "Peridiniella danica"              
## [153] "Amphidinium kesslitzii"           
## [154] "Gyrodinium spirale"               
## [155] "Prorocentrum cordatum"            
## [156] "Mesodinium rubrum"                
## [157] "Katodinium glaucum"               
## [158] "Chaetoceros debilis"              
## [159] "Attheya septentrionalis"          
## [160] "Pseudoscourfieldia marina"        
## [161] "Gymnodinium elongatum"            
## [162] "Thalassiosira nordenskioeldii"    
## [163] "Nematopsides vigilans"            
## [164] "Amphidinium crassum"              
## [165] "Apedinella radians"               
## [166] "Eutreptiella gymnastica"          
## [167] "Thalassiosira pacifica"
    # Plot
        multiOBIS <- sf::st_as_sf(multiOBIS,
                         coords = c("decimalLongitude", "decimalLatitude"),
                         crs="+proj=longlat +datum=WGS84",
                         remove=FALSE) %>%
            filter(!is.na(species))

        ggplot(multiOBIS) +
            geom_raster(data=bathyras,aes(x=x,y=y,fill=z))+
            geom_sf(data=bb,fill="transparent",colour='yellow',size=2)+
            geom_point(data=multiOBIS,aes(x=decimalLongitude,y=decimalLatitude))+
            coord_sf(xlim = c(lonmin-x,lonmax+x),
                     ylim = c(latmin-x,latmax+x),
                     expand = F)

Trophic network

    nsp <- length(spList)
    foodWeb <- matrix(ncol = nsp, nrow = nsp, data = 0, dimnames = list(spList, spList))

    pb <- txtProgressBar(min = 0,max = nsp, style = 3)
    for(i in 1:nsp) {
        for(j in 1:nsp) {
            inter <- rglobi::get_interactions_by_taxa(sourcetaxon = spList[i], targettaxon = spList[j], interactiontype = 'eats')
            if(nrow(inter) > 0) foodWeb[j,i] <- 1
        }
    setTxtProgressBar(pb, i)
    }
    close(pb)

    saveRDS(foodWeb, file = './foodWeb.rds')

We can now visualize the network using the package iGraph

    nsp <- length(spList)
    foodWeb <- readRDS('./foodWeb.rds')

    netw <-  igraph::graph_from_adjacency_matrix(t(foodWeb))

    vec_col <- spList
    vec_col[vec_col != sp] <- "#11bdec"
    vec_col[vec_col == sp] <- "#cc154b"

    plot(netw,
      vertex.color = vec_col,
      vertex.frame.color = "transparent",
      vertex.label.color = "black",
      edge.arrow.size = .5
      )

… Not very nice visually when there are too many species. We can however do this better (and more fun!) using the package networkD3

netwD3 <- networkD3::igraph_to_networkD3(netw, group = vec_col, what = 'both')
nodeSize <- spList
nodeSize[nodeSize != sp] <- 2
nodeSize[nodeSize == sp] <- 4
netwD3$nodes <- cbind(netwD3$nodes, nodeSize)

networkD3::forceNetwork(Links = netwD3$links,
                        Nodes = netwD3$nodes,
                        Source = 'source',
                        Target = 'target',
                        NodeID = 'name',
                        Group = 'group',
                        zoom = TRUE,
                        linkDistance = 50,
                        fontSize = 12,
                        opacity = 0.9,
                        charge = -20,
                        # colourScale = networkD3::JS('d3.scale.ordinal().range(vec_col)'),
                        Nodesize = 'nodeSize')

Neat huh?!?

Co-distribution

We could also evaluate how species are co-distributed using an analysis analog to SDMs, but at the communicty scale and using a hierarchical Bayesian approach, Hierarchical Modeling of Species Communities (HMSC; ref) currently being implemented in R in the package HMSC (ref). As with SDM using the biomod2 package, the first step is to format the data into a HMSC class object.

Formatting the data

Create regular grid over study area

In order to create common observations among individual observations, we will aggregate all data within regular grid cells, so that cells will become planning units on which we base co-distribution of species. We will also create a rather coarse sized grid in order to make processing speed more efficient. Of course there would be more sophisticated ways to do this, but for simplicity’s sake we will perform a simple intersect between the grid and occurrence data.

    # First create a function to generate hexagonal grids
        hexGrid <- function(size, shp_over) {
            library(rgeos)
            library(sp)
            shp_over2 <- gBuffer(shp_over, width = size) # buffer so that whole surface of shapefile is covered by a cell
            spat_grid <- spsample(shp_over2,type="hexagonal",cellsize=size, offset=c(0,0)) # creating points for a hex grid
            spat_grid <- HexPoints2SpatialPolygons(spat_grid) # creating polygons
            proj4string(spat_grid) <- proj4string(shp_over) # set projection
            spat_grid <- spat_grid[shp_over, ] # clipping grid to retain only cells intersecting shapefile

            id <- paste('ID',seq(1:length(spat_grid)), sep = '')
            row.names(spat_grid) <- id
            spat_grid <- SpatialPolygonsDataFrame(Sr = spat_grid, data = data.frame(ID = id, row.names = id))

            return(spat_grid)
        }

    # Load shapefile over which grid is made
        egsl <- rgdal::readOGR(dsn = "./", layer = "egsl")

    # Create grid
        egsl_grid <- hexGrid(size = 10000, shp_over = egsl)
        rgdal::writeOGR(egsl_grid, dsn = './', layer = 'egsl_grid', driver="ESRI Shapefile", overwrite_layer=TRUE)

    # Plot grid
        plot(egsl_grid)
    # Plot grid
        egsl_grid <- rgdal::readOGR(dsn = "./", layer = "egsl_grid")
## OGR data source with driver: ESRI Shapefile 
## Source: "./", layer: "egsl_grid"
## with 3172 features
## It has 1 fields
        proj <- "+proj=longlat +datum=WGS84"
        egsl_grid <- sp::spTransform(egsl_grid, CRSobj = CRS(proj))
        plot(egsl_grid)

Intersect occurrences and environmental data with grid
    # reload multiOBIS data and transform into sp points
        multiOBIS <- read.csv("multiOBIS.csv")
        proj <- "+proj=longlat +datum=WGS84"
        multiOBISpt <- sp::SpatialPointsDataFrame(coords = multiOBIS[, c("decimalLongitude", "decimalLatitude")], data = multiOBIS, proj4string = CRS(proj), match.ID = TRUE)

    # Overlay occurrences on grid
        occGrid <- sp::over(multiOBISpt, egsl_grid)
        multiOBIS <- cbind(multiOBIS, occGrid)
        multiOBIS_wide <- reshape2::dcast(multiOBIS, formula = ID ~ species, value.var = 'species', fun.aggregate = length)

        for(i in 2:ncol(multiOBIS_wide)) {
            for(j in 1:nrow(multiOBIS_wide)) {
                if(multiOBIS_wide[j,i] > 1) multiOBIS_wide[j,i] <- 1
            }
        }

        toAdd <- egsl_grid@data[,'ID'][!egsl_grid@data[,'ID'] %in% multiOBIS[,'ID']]
        matSup <- matrix(nrow = length(toAdd), ncol = ncol(multiOBIS_wide), dimnames = list(c(), colnames(multiOBIS_wide)), data = 0)
        matSup[,'ID'] <- toAdd
        multiOBIS_wide <- rbind(multiOBIS_wide, matSup)
        multiOBIS_wide <- multiOBIS_wide[-which(is.na(multiOBIS_wide[,'ID']) == T), ]
        egsl_grid@data <- merge(egsl_grid@data, multiOBIS_wide, by = 'ID')
        spAll <- colnames(egsl_grid@data)[2:ncol(egsl_grid@data)]
        spID <- seq(1:length(spAll))
        colnames(egsl_grid@data)[2:ncol(egsl_grid@data)] <- spID

    # Extract environmental values from rasters
        egsl_grid@data <- cbind(egsl_grid@data, raster::extract(bathypoly, egsl_grid, fun = mean, na.rm = T))
        colnames(egsl_grid@data)[ncol(egsl_grid@data)] <- 'bathy'
        egsl_grid@data <- cbind(egsl_grid@data, raster::extract(sst2, egsl_grid, fun = mean, na.rm = T))
        colnames(egsl_grid@data)[ncol(egsl_grid@data)] <- 'sst'

    # HMSC data
        studyGrid <- egsl_grid@data

        envCov <-  c('bathy','sst')

    # Remove environmental NAs
        NAs <- studyGrid[, envCov] %>%
                lapply(X = ., FUN = function(x) which(is.na(x))) %>%
                unlist(.) %>%
                unique(.)
        NAsID <- studyGrid[NAs,'ID']

        studyGrid <- studyGrid[-NAs, ]
Create HMSC class data
    # ---------------------------------
    # Y: sample units by species matrix
    # ---------------------------------
        # Extracting only presence/absence data
            Station <- studyGrid[, 'ID']
            Y <- studyGrid[, as.character(spID)]
            rownames(Y) <- Station

            for(i in 1:ncol(Y)) {
                Y[,i] <- as.numeric(as.character(Y[,i]))
            }

    # -----------------------------------------
    # Pi: sample units by random effects matrix
    # -----------------------------------------
        # Create a dataframe for random effects, columns have to be factors
            Pi <- data.frame(sampling_unit = as.factor(Station))

    # ----------------------------------------------------
    # X: sampling units by environmental covariates matrix
    # ----------------------------------------------------
        # Create a matrix for the values of environmental covariates at each sampling unit location
            X <- matrix(ncol = length(envCov),
                        nrow = nrow(studyGrid))
            X <- studyGrid[, envCov]
                    rownames(X) <- Station

    # ----------------------------
    # as.HMSCdata for HMSC package
    # ----------------------------
        # Creating HMSC dataset for analyses
            HMSCdata <- HMSC::as.HMSCdata(Y = Y, X = X, Random = Pi, interceptX = T, scaleX = T)
## [1] "column names were added to 'Random'"
## [1] "'Y' was converted to a matrix"
            str(HMSCdata)
## List of 3
##  $ Y     : num [1:2851, 1:167] 0 0 0 0 0 0 0 0 0 0 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2851] "ID1003" "ID1004" "ID1005" "ID1006" ...
##   .. ..$ : chr [1:167] "1" "2" "3" "4" ...
##  $ X     : num [1:2851, 1:3] 1 1 1 1 1 1 1 1 1 1 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:2851] "ID1003" "ID1004" "ID1005" "ID1006" ...
##   .. ..$ : chr [1:3] "Intercept" "bathy" "sst"
##  $ Random:'data.frame':  2851 obs. of  1 variable:
##   ..$ random1: Factor w/ 2851 levels "ID1003","ID1004",..: 1 2 3 4 5 6 7 8 9 10 ...
##  - attr(*, "class")= chr "HMSCdata"
HMSC analysis

We are now ready to perform the analyses!

    model <- HMSC::hmsc(HMSCdata,
                        family = "probit",
                        niter = 10000,
                        nburn = 100,
                        thin = 10)
    saveRDS(model, file = './HMSCmodel.rds')
HMSC analysis diagnostics
model <- readRDS('./HMSCmodel.rds')
# =========================================
# 3. Producing MCMC trace and density plots
# =========================================
        mixingMeansParamX <- coda::as.mcmc(model, parameters = "meansParamX")

    # Trace and density plots to visually diagnose mcmc chains
    # Another way to check for convergence is to use diagnostic tests such as Geweke's convergence diagnostic (geweke.diag function in coda) and the Gelman and Rubin's convergence diagnostic (gelman.diag function in coda).
        paramModel <- colnames(mixingMeansParamX)
        nParam <- length(paramModel)

        par(mfrow = c(nParam, 2), mar = rep(2, 4))
        for(i in 1:ncol(mixingMeansParamX)) {
          coda::traceplot(mixingMeansParamX[,i], col = "blue", main = paste('Trace of ', paramModel[i]))
          coda::densplot(mixingMeansParamX[,i], col = "orange", main = paste('Density of ', paramModel[i]))
        }

# ================================
# 4. Producing posterior summaries
# ================================
    # Average
        average <- apply(model$results$estimation$paramX, 1:2, mean)
    # 95% confidence intervals
        CI.025 <- apply(model$results$estimation$paramX, 1:2, quantile, probs = 0.025)
        CI.975 <- apply(model$results$estimation$paramX, 1:2, quantile, probs = 0.975)

    # Summary table
        paramXCITable <- cbind(unlist(as.data.frame(average)),
                             unlist(as.data.frame(CI.025)),
                             unlist(as.data.frame(CI.975)))
        colnames(paramXCITable) <- c("average", "lowerCI", "upperCI")
        rownames(paramXCITable) <- paste(rep(colnames(average), each = nrow(average)), "_", rep(rownames(average), ncol(average)), sep="")

    # Credible intervals
        paramXCITable_Full <- paramXCITable

        beg <- seq(1,nrow(paramXCITable_Full), by = nsp)
        end <- seq(nsp,nrow(paramXCITable_Full), by = nsp)
        sign <- abs(as.numeric(paramXCITable_Full[, 'lowerCI'] <= 0 & paramXCITable_Full[, 'upperCI'] >=0) - 3)

    # Export figure
        par(mfrow = c((length(beg)+1),1))
        for(i in 1:length(beg)) {
            paramXCITable <- paramXCITable_Full[beg[i]:end[i], ]
            cols <- sign[beg[i]:end[i]]
            par(mar=c(1,2,1,1))
            plot(0, 0, xlim = c(1, nrow(paramXCITable)), ylim = round(range(paramXCITable)), type = "n", xlab = "", ylab = "", main=paste(colnames(mixingMeansParamX)[i]), xaxt="n", bty = 'n')
            axis(1,1:nsp,las=2, labels = rep('',nsp))
            abline(h = 0,col = 'grey')
            arrows(x0 = 1:nrow(paramXCITable), x1 = 1:nrow(paramXCITable), y0 = paramXCITable[, 2], y1 = paramXCITable[, 3], code = 3, angle = 90, length = 0.05, col = cols)
            points(1:nrow(paramXCITable), paramXCITable[,1], pch = 15, cex = 1, col = cols)
        }
        mtext(text = spAll, side = 1, line = 1, outer = FALSE, at = 1:nsp, col = 1, las = 2, cex = 0.4)

# =======================
# 6. Association networks
# =======================
    # Extract all estimated associatin matrix
        assoMat <- HMSC::corRandomEff(model)

    # Average
        plotMean <- apply(assoMat[, , , 1], 1:2, mean)
        colnames(plotMean) <- rownames(plotMean) <- spAll

    # Build matrix of colours for chordDiagram
        plotDrawCol <- matrix(NA, nrow = nrow(plotMean), ncol = ncol(plotMean))
        plotDrawCol[which(plotMean > 0.4, arr.ind=TRUE)]<-"red"
        plotDrawCol[which(plotMean < -0.4, arr.ind=TRUE)]<-"blue"
    # Build matrix of "significance" for corrplot
        plotDraw <- plotDrawCol
        plotDraw[which(!is.na(plotDraw), arr.ind = TRUE)] <- 0
        plotDraw[which(is.na(plotDraw), arr.ind = TRUE)] <- 1
        plotDraw <- matrix(as.numeric(plotDraw), nrow = nrow(plotMean), ncol = ncol(plotMean))

    # plots
        # Matrix plot
        Colour <- colorRampPalette(c("red", "white", "blue"))(200)
        corrplot::corrplot(plotMean, method = "color", col = Colour, type = "lower", diag = FALSE, p.mat = plotDraw, tl.srt = 65, tl.cex = 0.4, pch.cex = 0.3, tl.col = 'black')

        # Chord diagram
        circlize::chordDiagram(plotMean, symmetric = TRUE, annotationTrack = 'grid', grid.col = "grey", col = plotDrawCol, preAllocateTracks = 1)
        circlize::circos.trackPlotRegion(track.index = 1, panel.fun = function(x, y) {
          xlim = circlize::get.cell.meta.data("xlim")
          ylim = circlize::get.cell.meta.data("ylim")
          sector.name = circlize::get.cell.meta.data("sector.index")
          circlize::circos.text(mean(xlim), ylim[1] + .1, sector.name, facing = "clockwise", niceFacing = TRUE, adj = c(0, 0.5), cex = 0.4)
        }, bg.border = NA)

    # ===============================================
    # 7. Computing the explanatory power of the model
    # ===============================================
        # Prevalence
            prevSp <- colSums(HMSCdata$Y)
        # Coefficient of multiple determination
            R2 <- HMSC::Rsquared(model, averageSp = FALSE)

        # Draw figure
            plot(prevSp, R2, xlab = "Prevalence", ylab = expression(R^2), cex = 0.8, pch=19, las=1, cex.lab = 1, main = 'Explanatory power of the model')
            abline(h = mean(R2, na.rm = T), col = "blue", lwd = 2)

    # =============================================
    # 8. Generating predictions for validation data
    # =============================================
        source('./crossValid.r')
        modelAUC <- crossValidation(data = HMSCdata,  nCV = 2, validPct = 0.2, as.strata = FALSE)
## [1] "The priors for the latent variables should be OK for probit models but not necessarily for other models, be careful"
## iteration 2000
## iteration 4000
## iteration 6000
## iteration 8000
## [1] "The priors for the latent variables should be OK for probit models but not necessarily for other models, be careful"
## iteration 2000
## iteration 4000
## iteration 6000
## iteration 8000
        meanAUC <- colMeans(modelAUC)
        sdAUC <- apply(modelAUC, MARGIN = 2, FUN = sd)

        plot(prevSp, meanAUC, ylim = c(0,1), xlab = "Prevalence", ylab = 'AUC', cex = 0.8, pch=19, las=1, cex.lab = 1, main = 'Monte Carlo cross-validation with AUC of ROC curves')
        abline(h = 0.5, col = 'grey')
        arrows(x0 = prevSp, x1 = prevSp, y0 = (meanAUC - sdAUC), y1 = (meanAUC + sdAUC), code = 3, angle = 90, length = 0.05)
## Warning in arrows(x0 = prevSp, x1 = prevSp, y0 = (meanAUC - sdAUC), y1 =
## (meanAUC + : zero-length arrow is of indeterminate angle and so skipped
        points(prevSp, meanAUC, pch = 19, cex = 0.8)

# Generate predictions
    HMSCpred <- predict(model)

# Replace NAs in grid
    data <- matrix(nrow = (length(NAs) + nrow(HMSCpred)), ncol = ncol(HMSCpred), data = 0)
    dimnames(data) = list(paste('ID',seq(1:nrow(data)), sep = ''), colnames(HMSCpred))

    for(i in 1:length(NAsID)) {
        data[which(rownames(data) == NAsID[i]), ] <- NA
    }

    for(i in 1:nrow(HMSCpred)) {
        data[which(rownames(data) == rownames(HMSCpred)[i]), ] <- HMSCpred[i, ]
    }

    HMSCpred <- data
    HMSCpred <- cbind(HMSCpred, rownames(HMSCpred))
    colnames(HMSCpred)[ncol(HMSCpred)] <- 'ID'
    # colnames(HMSCpred) <- spAll
    # Plot grid
        egsl_grid <- rgdal::readOGR(dsn = "./", layer = "egsl_grid")
## OGR data source with driver: ESRI Shapefile 
## Source: "./", layer: "egsl_grid"
## with 3172 features
## It has 1 fields
        proj <- "+proj=longlat +datum=WGS84"
        egsl_grid <- sp::spTransform(egsl_grid, CRSobj = CRS(proj))

        egsl_grid@data <- merge(egsl_grid@data, HMSCpred, by = 'ID')
        data <- egsl_grid@data[,100]
        rbPal <- colorRampPalette(c('#2f6eb9','#2aadba','#b45f5f')) # color palette
        cols <- rbPal(50)[as.numeric(cut(as.numeric(data), breaks = 50))]

        plot(egsl_grid, col = cols, border = cols)

<<BACK